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INTRODUCTION 


In this report, we summarize work done under NASA Contract 
NASl-15844 on compressible shear flows and drag reduction. The 
report consists of this Introduction and the three following 
sections that summarize our progress on different aspects of 
this work. 

In Section 1, we present a summary of some work on analytical 
and numerical aspects of conformal mapping. It describes a new 
and very efficient and robust method for computation of these maps 
in highly distorted geometries. 

In Section 2, we describe the computer code SPECFD written 
for the CYBER-203 at NASA Langley Research Center. This code 
solves the three-dimensional time-dependent compressible Navier- 
Stokes equations by a mixed finite-dif f erence-spectral , algorithm. 

It works efficiently on the CYBER-203 with resolutions of 32 x 32 x 
64 and promises to yield new insights into the nonlinear dynamics of 
compressible shear flows in wall-bounded geometries. 

In Section 3, we describe our v7ork on two-equation turbulence 
modelling of turbulent flow over wavy walls. A modified Jones- 
Launder model is used to include effects of the turbulence. This 
transport model has been implemented in a two-dimensional spectral 
code for flow in general wavy geometries. Results are presented 
for both flow over flat plates and flows over wavy walls. 


1 . 



1 . 


Analysis of Numerical Conformal Mapping 


Many two dimensional physical problems require the 
solution of Laplace's equation in a complicated domain . 

One way to solve these problems is to conform.ally map n onto 
the unit disk D(0,1) ^ plane D(~) • Once that is 

done the Poisson kernel provides the solution to the Dirichlet 
or restricted Neumann boundary value problems. Hilbert's 
■ generalization solves a mixture of the two where each applies 
on part of the boundary (but it doesn't solve the general 
Neumann boundary condition) . Conversely any method of 
computing the Dirichlet or Neumann solution can be used to ^ 
calculate the conformal map (see Theorem 5.3 ^of Dubiner 1981) 

but there is little reason to do it. 

There exists a unique conformal mapping f of n onto 

D(0,1) up to specifying f(v) and Arg 9v f (v) 

V en . Classical complex analysis demonstrates that on the 
boundary BJiyf is about as smooth as is and, of course, f 

is analytic inside. However, [8uf(u)] ^ is ill posed in 
terms of any reasonable norm of Q even when u is restricted 
to be well away from 3^ . For example take 


fi = - arctanh[tanh D(0,1)] 

IT ** 


(. 1 ) 


the notation means f2- {^arctanh [tanh ^ z]IzeD(0,l)) 


where 


2 


It is a smooth domain which looks like an ellipse inflated 
inside a rectangle centered at the origin of length 22, and 
width 2 - - arctanh(e'^) • But the conformal mapping taking 

TT 

n to. D (0,1) and 0 to 0 is 


f (u) = coth ^ tanh ^ 


so 


8uf (u) u=2, 
8^f (u)j u=0 


cosh 


--2 


1x2, 

4 


( 2 ) 


(3) 


which decreases exponentially in I and equals 0.000000603 
for t=10 '■ The curvature of on near the ends relative to^ 
n ’s diameter is OU) but it is innocent of (3); The 
eccentric cigar shape of S is to blame 'and the same would 

happen for the smooth paddle-shaped domain of Figure 1. 

Except near the ends example (.1) is a slender domain. 
A domain « is called e slender with c small (say, 0<e< 
iff 3n/{»} is composed of two connected components and 

such that for each 

1 I ~i ^ (4) 

|k(u,Wq) I . lu --al 1 e 


lArg 



I ,_< E 


(5) 


• 3 . 



where minimizes |u-u 1 <(.u^Wq) Wq ^ 

curvature at u. Condit ion ..(4) requires Wq to be nearly 
straight and condition ( 5) requires to be nearly 

parallel to Wq • Parametrize by its arc length s 

starting from an arbitrary fixed point. Each uefl can be 
■ uniquely written as 

u = Wq[s(u)] +t(u)n[s(u)] 0<t<h[s(u)] (6 


where n' (s) is the inside ^unit normal at WqCs) and hCs) is 
the distance ofn(s)from ^^(s) along the direction n(s) so 


that W(,(s)+h(s)n(s)eWj^ • Let us normalize 

system (s,t) in an approximately isotropic way 


the coordinate 
y 


' s(u) 

g(u) = 2 t / 

^ 0 


dp 

h{p) 


+ i 


t (u) 
h [s (u) ] 


The map g from onto A — {z] ] Im z| 
conformal with eccentricity bounded by 



C7)' 


< H} is quasi 

A __ ■ 

ce (see Dubiner 1981 


for definition)_^ g _be_the exact con fj^al map 

j sending Wq(±”) to ±®° respectively^ Clearly — 


f (u) = tanh g (u) 


(8) 


conformally maps 


Q onto D(0,1) and so does 



1-f (u)'f (v) 


(91 


and it sends ven to o . The real number c is determined by 


(v,v) > 0 


(10) 


. , « denotes differentiation with respect to the first 

variable. Formula (8) is inserted in (9) and results in 


f:, . _ sinh [g (u) -g (v)1 -i Arg3^g(v) 

f(U/V) cosh [g (u) -g (v) ] 


( 11 ) 


and 0-1) 's derivative is 


?ug(u) 

8 f(u,v) = 2 ' — 

^ cosh"^ [g (u) -g (v) ] 


cost2Im-g(v)Ie-^'''^'3Svg(v) 


Define f(u,v) by replacing 9 with 9 in (lO) . It is a 
quasi conformal map from onto D(0,1) of at most ce 
eccentricity sending v to 0 . Hence f(*,v) is expected to 
be close to f(*,v) in some sense. indeed (3.3) of Dubiner 
(1981) and others prove that 


i5,n 8u-5(u,v) - ilnSuf (u,v) I<c£ li(u)-g (v) I 


. (13) 


where 
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. Arg f(u,v) = - a [0 ,Wq [ s (u) ] ,W q] + 0(1) (16) 

where a[ • ] is the change in angle of between 0 =Wq[s(v)] 
and u»s projection on , Wq[s(u).] , Thus globaly f 

performs reasonable rotation but extreme scaling. In 
retrospect it should not be surprising because conformal maps 
are defined by being locally angle preserving with no scaling 
restrictions attached. 

Formula (l6) is easy to interpret. It obviously holds 
(up to translation in Arg 3^f (u) depending on its 
normalization) for ue3f2 r where n is a general domain. Thus ' 
(16) states that for slender domains 

I ^ 

Arg 3^f(-,v) ^ = 0(1) (17) 


7 . 



where the notation means ]Arg 3~f(Q,v) - Arg3^f(u,v) l£c and 
ue3n is near u , say the closest boundary point. The result 
(17) holds in general as proven by Theorem 5.4- of Dubiner 
Ti981)'. ' Formula* (15) is not that easy to generalize. Unlike 
(16) , its right side depends on the structure of ^ betv;een v 

and u The f ir s t gue s t ion _i s j wh at does 'between* m ea n in 

general? In order to gain some insight let us consider a more 
complicated example. 

Let 0 < .£ << 1 

fi(e) = (x+iy 1 ex + cos y > 0 } 

The domain (18) has the following property. . Any domain ^ 
is said to be a c > 0 conjugation of the domains' iff 
for any uefi there exists a vel and two complex numbers a,b 
such that 

ueA (u) = aA^ + b (19) 

a[u,A(u) £ e (20) 

the distance from A(u)to n relative to u is defined by 
(10.113) in Dubiner (1981). The interested reader may prove 


that any e slender domain is a ce fconjugation of 
. = {x+iy] Jy] < J-} 


( 21 ) 


where oO is constant. Domain (18) is a ce conjugation of 
and 

A2={x+iy|x - 0} ^ (22) 

A^= c\{-A 2 ) = {x+iylx+ > 0 } (23) 

CO 

A- = c\u [-“,i (2n+l) it] (24) 

n=-oo . • 

where [a,b] is the closed interval betv;een a and b . V7e 
have to match the conformal maps from all the A(u),ue^(£) 

In this case it is easiest to do when considering f [ * (e) ], 

the conformal map from fi(e) onto the half plane D(«>) 
normalized by 

8^f (■+“,“) =1 (25) 

The domain J2(e) is periodic and symmetric so we can limit 
ourselves to 


9 . 



( 


u = x + iy 0<.y<7T 


(25) 


We start from 


■u 


f(u,oo, A^) = 2 arcsinh e' 


(27) 


It is modified to 


X > — - 0 (— ) |x - — + i (y-ir) I >> e 


(28) 


f[x + iy / J2( e) ] = 2 arcsinh e 


1 , ,1. . . 77 

5(x - -)+ ly 


2h (x) 


(29) 


where 


h(x) = \ axccos {—cx) 


- 1 < =. < i 

e — — e 


(30) 




!<>= 


(31) 


and the exact limitation (28) will follow from comparison 
with the following formulas. A priori rigorous bounds can be 
derived but as usual it is inconvenient. Next 


f(u,°°,A^) = i c e 


2u 


(32) 
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(there is no natural normalization)*. V7e already know hov/ to 
match the f[u,”,A(u)] . ’s where A(u)=aA^+b : recall 

(7,8). They match with (28,29) and give 


(33) 


|-< X < ^ + 0 (— ) V 


[x + i + iyl>>£, lx- |+i(y-iT) !>>e 


f [x+iy (e) ] = 2e 


irfsi rh (x) ]-si (tt) . y I 

2L i h(x)J 


(3A) 


where 


SI 


(p) = I 


sin q 


dq 


(35) 


Now we can match f(*,“,A^) to (33,34). and obtain 


jx + i + iy 1 


-1/3 


<< e 


(36) 


TTSi(lT) 


2e 


cos 


- |[x + I + iy 


f [x+iy,“,n (e) ] = 4e 
Similarly, f(-,”,A 3 ) is matched to (28,29): 


/2EX+2 , 

h(x) 1 
(37) 


|x + i + i(y-TT) 1 < < 1 


(38) 


12 


f [x+iy (e) 1 = /2z ^ [x - ^ + i[y-h(x)] ~ 

( 39 ) 

In particular the maximum and minimum of I 9 ^f [u,<»,Q (e) ] [ 
are obtained at ^ + iir , - respectively and 

^ + iTT,co,n (e) ] ^ /|T (40) 


TTSi(lT ) 

9lf [ - I ,»,n(e)]- ^ e 


(41) 


What have we learned from example (18)? Figure 3 
illustrates the direction of information flow (The reverse of 
the direction of dependence) which were exhibited while 
fl»,“,n(e)] has been constructed. The situation is quite 
special yet v;e have some grounds to suspect that in general 
9^f(u,v,f2) and other functions depend mainly on ' s part 
’around' the curve of least Euclidean distance betvreen v and 
u inside fi. A close inspection of (27-38) reveals that 
the above mentioned curve from to u resembles 


r(~,u,fi) 


{wef2 I 


n f (w,”,n) 
• f(u,“,fi) 


< 1 } 


(42) 


The curvces r(u,v,fi) are called geodesics because they are 
the geodesics of a certain ccnformally invar iont metric 
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far 


p(u,v,n) of Theorem 1.1, Some geodesics of are 

illustrated in Figure 0.4. Notice for any tv;o u,veA^ 
away the most of (u,v,Aj^) is exponentially close to A^' s 
axis of symmetry. Theorem 8.3 demonstrates that in general 
the geodesics try to keep away from the boundary. 

The connection between geodesics and lines of least 
Euclidian distance is proven in Theorem 9.2 (Dubiner 1981). 

Now that we have some idea on what 'between’ and 
means it is time to find out how the rest of n affects 
An 8^f(u,v,f2) and other quantities. For that purpose 

let us return to e slender domains and compute a next order 
correction to (7) . We start by calculating 8^ . The 
gradient of (6) is 

du = n [s (u) ] [ (-i+t (u) 3 An n [s (u) ] ) ds + dt] 

s 

and of (7) 

dg(u) = 2 hW [ (l-it{u)3gAn h(s))ds + idt] 


(43) 


( 44 ) 


Formulas 


(43,44) combine into 


7T 


dg 


4hn 


•t3 iln(hn)* 

[2i + — -] 

l+it3 J-nn 
s 


du - 


n 

4h 


t3 S,Ti (hn) 

— ^ du 

l+it3 £nn 
s 


Define the g correction function 


( 45 ) 


g^-= g - g 

Then 


8_'g, = - 3_ g = - 2iri & y= ” 27 rir| 

u u 


Im g^l an = 0 


where - 

. , it(u)3 S,n(h[s(u)]n[s(u)]) 

y(u) = I ^ 

1+it (u) 9 g 2 ,nn [s (u) ] 


Notice that is 0(e) slender iff 


(46) 


(47) 


(4B) 


(49) 


sup I y (u) I ^ ce 

uen 


( 50 ) 


Problem (47,48) has a unique solution .up to an additive real 
constant ■ 

. .1 .. I — — — - , % 

(51) 


. ^ [■ n (w) 3^f (w) n(w)3^f(w)-| 

” ^ 1 f (u)-f (w) 1 _ 

“ fin) 


g^(u) = 


f (W; 


I I 


g(u) = g( 


u) = i / J n (w) 3^ q (w) coth [g (u) -g (w) ] - 

i'2 L ^ 




where 


n (w) 3 g (w) 

W .r 


tanh [g (u) -g (w) 


Id^w] 


(52) 


|d^wl = d Re w* d Im w 




(53) 


Formula (51) is an integral equation of the first kind which 
can be iterated to convergence. The first order correction to 
g with some modifications is 


g(u) ~U + / / [y (w)coth (U-W)+TTTwO tanh(U-W) ] Id w| 

n 


(54) 


where we have and will abbreviate 


U = g(u) , V = g(v) , w = g(w) 


(55) 


Recall (11) 
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( 56 ) 


Itn t(u,v) = In 1 to ^ 

cosh (U-V) ^ d V 

so 


Jtn f(u,v)^ In f(u,v) + • 

+ / / lu(w)K(U,V,W)+ir(^)K(U,V,W+ij) ] |d^wl 

n 

K(U,V,W) = coth(U-V) [coth(U-V7)-coth(V-T'Jl] - 
- tanh (U-V) [coth (U-W) -tan 1 (V-W) ] + 


+ i sinh“^ (V-W) - ~ sinh"^ (V-W) = 


' cos (2lm v) p Sinh (U-V) 

2sinh (U-W) sinh (V-W) sinh (V-W) ^sinh (U-W) 


and 


£n 8^f(u,v) = £n f(u,v) + £n 3^ £n f(u,v) 


so 


(57) 


cosh (U-V) 
cosh (V-Vv) 


(58) 


. (59) 
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■ (61) 

If we assume that U(w) is not only of order e but is 

also slowly varying, for instance — I 9 .U [w„ (s) +h (s) n (s) ] ! 

. n (s ) s 0 — 

then the integration in (57,60) can be done explicitly. Of 
course that results is much easier to derive directly and is 
of no interest to us. What we have wanted and obtained is the 
relative dependence of 3^f(u,v) on the boundary part 

a= Wq([s- iAs,s+ |a s]) centered at w = WQ(s)of length as . 

When 

0 < AS £h(s) £ |u-w| , |v^w] (62) 

the contrlling factor of “ ' s influence is 

g-4infig(w) - (g(v),g(u)!!3up|„ (a) 1 (63) ^ 
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The asymptotically correct term is the same with g replaced 

by g. A graphic interpretation of (63) follov;s. In order 

to affect 3 f(u,'v) the data about a's shpae must travel 
u 

from a to V . It is provided a free ride in Q 's portion 
betv/een u and v (in general on the geodesic F (u,v) betv/een 
u and v ) but it must pay 21iTs (zTT distance for travel 

around any point z not between u and v . The data is 
thrifty so it will move along a geodesic (a cost minimizing 
curve which turns out to agree with ( 42) to some point 0 
between v and u and will then enjoy a free ride to v . The 
optimal choice of 0 is the middle point among u,v,w . Let 


us call the total minimal cost p (w,As,u,v) 


Then v;hen 


the a data reaches v its intensity is deminished by a factor 
-4Px 

of e . A similar situation holds for other function of 


the conformal map besides 2.n3^f(u,v) , except that the ride 

on the geodesic betv?een u and v is not free but on a reduced 
fare. For instance (57,58) shows that the controlling 
factor of a's influence on in f(u,v) is 


^-2inf I5(w)-[g(v) , g (u ) ] ! -2 1 g (w) ~g (v) Ig | 

which means that the travel on F (u,v) is done on a half fare. 
Notice that this, does not change 0 , the point of transfer to 
f(u,v) 

Nov; that we know v;hat ( 63) means let us understand ' 
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-2P 

where it comes from. The term e ^ i& simply the decay of a 

Dirichlet or Neumann data from ct to f(u,v) and is already 

present in (54) ; Hov/ever, in the computation of (58) a 

-2P 

X 

cancellation has occurred and another e has appeared. 

Lest it look like a freak accident let us derive it from 
another point of view, close in spirit if not in technical 
complexity to section 10 *s. Suppose that 0 , is not only e 
slender but is e close to A. If w is between u and v 
than p = 0 so let us consider v between u and w (the 

X 

remaining case is similar) . Define 

i f(z) = tanh(z ~ Re v) (65) 

I ■ ' ■ ■ 

i ■ ■ ' ^ ^ 

i The domain f(S7) is close to the unit disk. The image of a 

i- "2Px 

I has length of order e As . A major change in a modifies 

i f(f)) by a region av-zay from ? (u) , ^ (v) whose area is of 

j order (e ^As)"^ . The details will be presented in Section 

10 but it takes no great leap of imagination to conclude that 

1 

i a 's effect is at most proportional to that area, and that is 

! -2p 

j where the extra e ^ comes from. What about g(u) ? It is 

; normalized so that g(±«>) = 1 and ] ? (w) - Ij or | f (w) + 1] 

j ' -2p 

is of order e ^ so '^'s influences on g(u) is of order 

~2P„ 9 

(e ^As)^ 

q“2Px 

The interpretation of (63,64) v/as chosen so that it it 
: general izable to arbitrary domains, v/ith some modifications. 
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The asymptotic theory has to be replaced by estimation up to a 
constant factor. The measure V of. ^'s deviation from A «s 
shape is special and anyv;ay there is no ideal general domain. 
Instead v/e v/ill pick a wide challenge and consider a 
perturbation of a general domain fi to a specified general 
domain n . We will divide the perturbation into parts and ^ 

as in Theorem 10.5 I show that each such part a of diameter As 
centered at w af fects S.n a^f(u,v) by 6^(w,AS,u,v) which 

is a generalization of (63) with sup|y (a) | replaced by 
The term 5(w,As,'u,v) is first encountered in 

Theorem 8 . 6 and 

6 (w,As,u,v) = sup 6(w;As,u,v) (66) 

^ der(u,v) 

The general interaction betv/een u,v,w and 0 is described 
by Theorem 4.6, which should be combined v/ith Theorem 9.2 

(Dubiner 1981) . 

The local length scale h[s(w)] is generalized into 
a(w,u,v) of (9.2). Its dependence on u,v is unfortunatly 
unavoidable: consider a half planed D(«>) . The lack of a 

local length scale, except of inf 1 w-9n| which vanishes on 
the boundary is the most imtportant general ingredient 
missing in slender domains. It is not fortuitous that the 
conformal m.otri’c. p(u,v) blcws up at the boundary. 

Our examples were mostly of smooth domains but notice 
that domain ( 18) has a vc; • arp bend at ^+i"^ , v/hich 
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rresponds to A 's corner, and*it did not distrub us from 
completing a uniform asymptotic approximation. Theorem (10.7) 
proves an interesting property of a fractal, which is the 
applied mathematician's ultimate in roughness. Ho\;ever, 

~we insist' on obtaining specific estimates at specific points, 
in contrast to the 'average' type approach of P.D.E. Theory. 

That is an advantage when it works but it fail^ hear a jrough 

boundary. In Sec. 5 of Dubiner (1981) it is shown hov; to 
patch our results with P.D.E. Theory. Some synthesis is clearly 
required. 

Let us now consider the numerical computation of 

conformal mappings. As remarked earlier, any Laplace solver 

• 2 

will do. Suppose that the domain is covered by 0 (N ) 

points, N of which are on the boundary . Then a Laplace 

2 

solver requires storing 0 (N ) numbers and performing from 

o 3 

0(N .S-n N) to 0 (N ) operations, where the last estimate is 
more realistic for complicated domains. The grid set up is 
is troublesome, especialy for multiscaled and time dependent 
(free boundary) domains. 

One way to avoide an internal grid is using a vortex 
representation. That results in an integral equation of the 
first kind which is numerically formulated as a set of n ^ N 
linear, equations. It can be solved by by Gaussian elimination 
which requires 0 (N ) memory locations and 0 (N ) operation. 
Alternativly one may iterate the system using 0 (n) memory and 


0(N ) operations per iteration. This is the numerical approximation 

to Neumann's series, and the later is guaranteed to converge for any 

single- sheeted domain satisfying some mild conditio.ns. The best 

existing Rayleigh-Taylor, instability simulation has been done in that 
way by Baker, Heiron and Qr saag (1931) . ■■ It is relatively easy to program 

and generalize to 3 dimensions. Moreover it can handle two incompressible 

fluids problems which conformal mapping alone can not solve unless the 


fluids density ratio is O' or 1. Hov/ever, Neumann's series's convergsnci 
is precarious. The arte of convergence for the domain ^ equals that for 
its exterior domain c Q and in particular convergence fails for multi- 
sheeted domains and is very slow v/hen two separate parts of approach 


each other. The domain 


(18) 


requires 0(g) iterations per order 
reduction in the error. Moreover thb Neumann series seems hard to modify 
in a way which v;ill extract a singularity such as a corner and still preserv 
convergence for general domains. 

The most natural numerical conformal mapping computation is done 
by Taylor expanding the conformal function from the unit disk onto ^ . 

Several such methods are listed in Sec. I. The best of them takes only 


0(N) ; memory locations and 0(N^£n N) operations but we have seen that the 


series will not converge to domain (is) 


before N = e 


c/fe 


terms are taken 


The first direct computation of the conformal map onto a 
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cigar shaped domain has been done by Manikoff and Zemach (1980) . Their 
method is to set up a system of N nonlinear equations and .solve then b 
Nev/ton iteration. Each iteration takes O(M^) memory and O(N^) 
operations and only fev; iterations are required. 

Any partial different ial equ ati on on a time dependent domain 
can be solved by a Green’s function method which utilizes only 

2d-2 

boundary data. This is not usually done because it takes 0(N ) 

memory locations and operations in d dimensions which is 

unreasonable for d ^3. For the operator in 2 dimensions better can be 
done because then the Green's function G(u,V/^) is constructable from 

from G(u,VQ,fi) and its harmonic conjugate where is 

constant. ' That is the basis of our method, though it will be 
presented in a different way. It ta)Ces 0(N) memory locations 

and 0 (N^) operations. 

1 \ 
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2. NUMERICAL SIMULATION OF COIIPRESSIBLE SHEAR FLOWS 
ON THE CYBER-203 


INTRODUCTION 


The objective of this work is the .development of 
an algorithm for the numerical simulation of three- 
dimensional viscous compressible flows. This code 
is applied to basic studies of compressible shear 
flows at high Reynolds number, such as the simulation 
of the incipient stages of transition to turbulence, 
and receptivity of laminar' boundary layers to 
external disturbances. Since the physical- instab- 
ilities involved in such problems are of a delicate 
nature which could be easily masked by numerical 
instabilities or other errors, it is essential to 
use an accurate numerical algorithm that can be 
efficiently implemented on the Control Data Corp. 
CYBER-203. 

The primary purpose of this paper is to describe 
the mixed spectral-finite difference code, SPECFD, 
as coded for the CYBER-203. We will focus on the 
vectorization techniques for some of the computational 
procedures with careful attention to the interplay 
between the storage allocation, the vectorization 
techniques, and their combined effects on input/ 
output (I/O) requirements for this virtual memory 


machine. 


The next section describes the numerical method. 

We then discuss the characteristics of the CYBER- 203 
at NASA Langley Research Center. The next section 
describes vectorization techniques including data 
base design and computational procedures. Finally, 
the overall flow of the code is outlined and a timing 
analysis of the various parts is presented. 


Governing Equations and Solution Technique 


The compressible Navier-Stokes equations are 



= - q.Vp-pV*q 


Bp _ 


q* Vp + r (p 


av 





Pr 


yv 



+ 


(y-l)$ , 


where q, p,p,x,$,a,Pr is the velocity vector, density, 
pressure, viscous stress tensor, viscous dissipation, 
ratio of specific heats, and Prandlt number, respectively. 
The constitutive relation for the stress tensor is 


T = X V • q I + y (V q + V q ) 

2 

where X is assumed equal to - -jy and y is 
given by Sutherland's formula 

y = 1.4582 x'lO“^ T^^ V (T +.110.33) p^-sec 
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' The expression for viscous dissipation is 

$ = X(V-q)^ + ^ [(Vq) + (Vq)'^] * [ (Vq) + (Vq)'^]. 

The^ assumption of a calorically perfect gas is 
implicit in this formulation. 

The addition and subtraction of the term (1/p) Vp 

dv 

in the momentum equation and a similar manipulation 
in the equation for pressure require explanation. 

(1/p)^^ is the average of the specific volume on 
each horizontal x-y plane and hence is a function 
of the z-coordinate only. The Crank-Nicolson method 
is used on the pressure term in the curly bracket. 

We resort to this kind of artifact to ensure convective 
stability at low Mach numbers without having the severe 
time step restriction demanded by an explicit method. 

Boundary Conditions 

We choose periodic boundary conditions in the 
streamwise (x) and spanwise (y) directions. This 
permits spectral representation of the dependent 
variables in these two directions. In the normal (z) 
direction we use no-slip boundary condition on the 


plate (z = 0) . Using C = az/(az+b), the serai- 
infinite physical doraain 0 _< z < is raapped onto 

the coraputational doraain 0 £?< 1. Zero perturbation 
boundary conditions are applied at ^ = 1. For free 
shear flows^ the physical doraain - co < z < is 
mapped onto 0 £ C <_ 1 using a hyperbolic tangent 
transformation. Zero perturbation boundary conditions 
are applied at both top and bottom boundaries . 

Numerical Method 

The numerical method consists of three stages 
(fractional steps) based on the well-known technique 
of operator splitting [1]. The first fractional 
step includes the effect of the advection terms: • 
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p = + (Atj^ + At^) [-q^'Ap^ - p^V*q^] 


r -^n-1 „ n-1 „n-l„ -tn-l, 

- At^[-q *V ,p - P V*q ] 


p" = p'^ + (At^ + At'')[-q''-V p"" - r(p'"-p^^)V-q'"] 


- At^I-q *Vp - r (p - -P^-y )V’q ] 


— — At^ p V*q ^ 
‘ 2 ^av ^ 


where At^ = (At^)^/(2At^ ^) . A variable-step 
second-order Adams-Bashforth method has been used for 
time discretization. The pseudospectral method is 
used for calculating x- and y-derivatives , and 
central differences are used to evaluate z-derivatives 
(using an equally spaced grid in the 5 variable) . 

The second fractional step is an implicit 
pressure correction; 



This system of equations is reduced to a single 

* * • 

Helmholtz equation for p : 


** 

P 



Pav’-"4> Vp**J = P* 
p av 



r p V 
av 


(3) 

* * 

Using a spectral representation for p in the 
X— and y-directions and central finite difference 
in the z-direction there results a tridiagonal matrix 

‘ "k "k 

equation for p in wave number space. 

In the third fractional step / viscous corrections 
are added to the inviscid solution obtained in the 
previous two steps. Use of the pseudospectral technique 
for the evaluation of viscous terms entails a rather 
unreasonable increase in the number of FFT's per time 
step [2]. Consequently, the viscous terms are 
discretized using central differences. The truncation 
error is on the order of v ( A t + Ay + Ay + Az ) 
where v is the maximum kinematic viscosity. This 
error is small since v is small in the high Reynolds 
number flows of interest. Further, the well-known 
problem of artificial viscosity (due to t.he truncation 


* 

•q = Q 
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error of the advection terms) swamping the physical 
viscosity does not arise in the present algorithm, 

V 

as the advection terms are evaluated pseudospectrally 
To avoid severe time step restrictions resulting 
from small Az in the stretched z-mesh, the z- 
derivatives are treated implicitly. This results 
in tridiagonal matrices for u,v,w, and p. 

We write this three level time split scheme 
symbolically as 

U^+l = ^ 

where U = (u,v,w, ,p)"^ and T indicates transpose. 
Here the subscript E indicates the explicit 
step, I the implicit pressure step, and V the 
viscous step. 



Characteristics of the CYBER-203 


The CYBER-203 is a vector processing computer 
manufactured by the Control Data Corporation (CDC) . 

As a vector processor, it can perform vector operations 
on 64-bit (32 bit) data at a peak rate of 50 (100) 

MFLOPS (millions of floating operations per second) 

25 (100) MFLOPS for multiplication, and 12.5 (50) 

MFLOPS for division. With each vector operation 
there is a vector startup which is independent of 
the vector length. Vector operations of length 
1000 achieve 80 to 85% of their peak rate. This 
is about four times faster than vectors of length 
75. Consequently, it is desirable to use long vectors 
where possible. The ability to do this depends upon 
the algorithm being used as well as the storage 
scheme selected by the program designer. The CYBER-203 
can also do scalar processing and in this mode generally 
will execute two to three times faster than the CDC 
CYBER-175. The CYBER-203 has a virtual memory 
implemented as "small" pages (512 words) or "large 
pages (65,536 words). The CYBER-203 at Langley 
Research Center has 1025K, 64-bit words of semi- 
conductor memory (16 large pages) . One can effectively 
double the memory capacity by using the 32-bit word 
length features. 
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Vectori 2 ation Techniques 


The data base design and vectorization techniques 
are inter-related on the CYBER-203 since vectors, by 
definition, are contiguous locations in memory. 

Notation: 


We first introduce some notation helpful in 
describing these two facets of the implementation. 
Assume that the dimensions of the grid ‘in the 

x,y,z: directions are I,J,K. U will refer to the 

' T 1-9 

vector (u,v,w,p,^) and an operator such as 

operating on D should be applied to each component 

of U.. References to a "plane" of data are to one 

of the x-y planes. These are typically the vectors 

used in our code. U(k) refers to the kth plane. 

The Thomas algorithm for solving a tridiagonal system 


requires the storage of two pieces of data for each 
equation, and consequently each grid point: (a^(k),B^ 

are the kth plane of such data for the w-momentum 
equation. Also, U* is the result of the 
operator on l/', U** the result of the ^ 

operator, and the result of the operator. 


(k)) 


Data Base Design; 


We wish to accomplish two things with our data 
base design; (1) a storage scheme that permits operation 
on long vectors^ preferably over an entire plane of 
the three dimensional grid, and (2) a scheme that 
allows controlled paging by the virtual system. With 
this in mind, each variable which is required over the 
entire grid is stored on one x-y plane at a time, 
with all variables at the kth plane stored in 
succession (see Fig. 1) . The execution flow is then 
to do all computations possible at plane k before 
proceeding to plane k+1. It has been shown [3] 
that this data layout and computational strategy is 
important for controlling paging when the large pages 
are involved. The paging per time step is further 
reduced by (1) marching forward, then in reverse 
through the data base (z-direction) and (2) applying 
all three operators in one forward-reverse sweep. 

The first technique permits the I/O requirements 
for each time step to be proportional to the amount 
by which the data base exceeds central memory [3] 
and the second technique means that only one pass 
through the data base is required per time step 
instead of three. The latter approach does cause some 



programming inconvenience when computing the z- 
derivatives. For instance, the operator computation 

i( * ic 

at the kth plane requires u (k+1) - u (k-1) , 
where u is computed by the operator. The 

£i 

* 

requirement for u (k+lj implies that the 
operator must be lagged one plane behind the 
operator. 

Main Computation Modules; 

There are three algorithms which are a major part 
of the overall computational procedure and which have 
been vectorized for the CYBER-203. These are: 


1. One Dimensional Fast Fourier Transform (1-D FFT ) 
The X and y derivatives in the and 

g operators are computed using the one-dimensional 
FFT. For instance 




^ u(x). = F ^ (iug(oj) ) , 


where g(w) is F (u) , the Fourier transform of u 
with respect to the x variable and m is the wave 
number. The CYBER-203 implementation described in 
reference [4] is used and will be detailed more 
fully under the section "Performance". 
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Two-Dimensional Fast Fourier Transform (2-D FFT) 


The operator involves applying the two- 

dimensional FFT operator, to equation (3) . 

If we define 

P** (k^,ky,z) = F^yP** (x,y,z) , 
then the resulting second order differential equation 

"k k 

in z for P can be central differenced leaving 

an independent tridiagonal system, implicit in the 

z-direction, for each (k ,k ). wave number. The 

X y 

right-hand side for each system is the appropriate 
element from ^^y^* ' coefficients can be 

analytically determined. The solution to each system 

k k 

yields P (k^,ky,z) and then one computes- 
P**(x,y,z) = F"y P(k^,ky,z). ■ 

The CYBER-203 software described in reference (5) 

is used to apply F and F ^ . 

xy xy 

3 . Systems of Tridaigonal Equations 

The and operators both generate systems 

of tridiagonal equations implicit in the z-direction 


There is an independent system for each (k^,ky) or 
(X/Y) point, respectively, in a plane. The usual 

m 

Gauss elimination based tridiagonal solver, frequently 
called the Thomas _alspxithmt^is inherently scalar but 
since we must solve, say , M - IxJ systems 
simultaneously, vectors of length M can be used. 

The forward sweep of the Thomas algorithm is done 
a plane at a time saving only the pairs (a (k) , 3 (k) ) 
required for the back substitution of a unit upper 
bidiagonal system of equations. Since the variables 
are interleaved a plane at a time, the storage is 
appropriate for both the Iqng vectors and the 
controlled paging. 



Performance 


The CPU efficiency of any highly vectorized program on a 
vector computer like the CYBER 203 relative a pa'rticular scalar 
computer is directly proportional to the vector lengths involved 
This is true assuming: 

1. The same algorithms are implemented in both the scalar 
and vector computers. 

2. No extra work, typically in form of data movement, is 
required to achieve the vector lengths used. 

Program SPECFD is highly vectorized. For an I x J x K 
grid all vector operations, except the evaluation of the one 
and two dim.ensional FFTs are of length I*J or I*J+2I and, hence, 
quite efficient. The derivative'evaluation by the FFT routine 
is somewhat more difficult to analyze. 

The software to implement the FFT algorithm simultaneously 
computes the one-dimensional FFT of M independent complex data 
sets each of size N. It requires that the first components of 
each data set to be stored consecutively, followed by the second 
components, etc. With this arrangement, it achieves average 
vector lengths of Mlog N. It is obvious that the data structure 
is quite important and that increases in M will yield higher CPU 
efficiency. 



iLe total data base includes 17 variables stored over the 
entire grid. These are (u,v,W/p,p’) , their respective convective 
derivatives (durdvrdwxdpfcip ) from the previous timestep for 
the explicit operator, and the four Thomas algorithm pairs 
(ttu r3u )r (“v '^v Since 

only one of these is stored. For a 32 x 32 x 64 grid , the 
data base, code and ^stem requirements are nearly 21 large 
pages. Testing has shown that each forward-reverse sweep 
gives only 5 large-page faults per time step, as predicted. 
SPECFD uses the CYBEia. 203 ADVISE feature to overlap the I/O 
time with CPU execution. The total CPU time for a 
run is about 3.1 sec per time step. 


32 X 32 X 64 


3. Niimerical Simulatibn of Turbulent' Flow Over Wavy Walls 


Introduction 

A general computer code to solve the Navier-Stokes 
equations in two-dimensional geometries with a 
turbulence model has been developed. Here we report 
on the results for flows over wavy walls including 
turbulence effects by means of one- and two-equation 
turbulence models and compare the results with 
available experimental results. 

In order to compute turbulent flows in complicated 
geometries at realistic Reynolds numbers, it is necessary 
to model the effects of turbulence. Perhaps the 
simplest model is a one-equation model in which the 
Reynolds stress is determined by the mean-velocity 
gradient in terms of a relation of the form 


-(u.u. --=-5..k) = 
3 i: ' 


au. 

3 


au. 

1 


where is an appropriate eddy viscosity. The 

eddy viscosity is then modeled in terms of a 

turbulent length scale L as 


V 

T 



U 

a X 


au. 

'ax. ax. 

3 1 


9u^ 1/2 

)] 
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The turbulence length scale 'L' includes effects 
of pressure gradients as well as other effects such 
as the presence of a wall. The choice of ’L' 
will generally involve one or more free constants 

that must be chosen to match experimental results. 

One equation turbulence models of the form 
(I.l) have been quite successful at modeling flows 
that are in equilibrium. However, the constants 

used in the closure must be changed when the flow 
changes. In highly non-equilibirum flows (such as 
flows undergoing strong acceleration) , the optimal 
choice of these modeling constants may be difficult, 
if not impossible, to achieve. Also, in flows where 
there is flov; reversal, the utility of one-equation 
models is dubious. 

In the last decade, there have been several attempts 
to improve upon the one-equation models of turbulence 
by means of multi-equation models. In two-equation 
models of k-e type, gradient transport equations are 
developed for turbulent kinetic energy (k) and 
turbulent energy disipation ( e ) . Such two-equation 
models have been successfully applied to the calculation 
of flows in the presence of strong pressure gradients. 

There has also been some work with three- and four-equation 
turbulence models but the payoff for the increased computational 
and analytical complexity of these latter models has not 
yet been spectacular. 



The principle drawback to the use of two-equation 
turbulence models is that it increases the computational 
complexity of the system of equations one has to solve. 
There is also the difficulty of having more free constants 
in the model than in one-equation models. It is only 
in the past few years that these two-equation models 
have been coupled with full two-dimensional Navier- 
Stokes codes. The present work is notable in that 
we present results for the time dependent Navier- 
Stokes equations coupled with two-equation models in 
complicated wavy geometries. We choose to use a 
modified version of the two-equation turbulence model 
due to Jones & Launder (1972), as modified by Chien 
(1980) . The test results discussed below are presented 
for relatively low Reynolds number flows because of 
the special interests of these flows in applications. 

At low Reynolds numbers, the validity of the two-equation 
models is severely tested. 
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Method of Solution 


The two dimensional Navier-Stokes equations are 


8v 

a 

9t 


+ 


V Vv 

a 


= _ 1p_ 

9x 


+ 


a 


9 



(II. 1) 


and 


V*v = 0 (II. 2) 

where the stress tensor T includes both viscous 
stresses and turbulent Reynolds stresses. Equations 
(II. 1, 2) are solved in the region 

0 < X < 2tt, f (x,t) £ y < 00 (II. 3) 

above the wall y = f(x,t). Periodic inflow-outflow 
boundary conditions are applied in x: 

v(x ;< 2it, y,t) = v(x,y,t) (II. 4) 

A conformal mapping technique (Meiron,Orszag & 
Israeli 1981 and Sec. I) is used to transform the 
region in (II. 3) into the region 


0 < Z < 2it, 


0 < n < 


00 


(II. 5) 


non-conformal stretching of the n variable is 
used to implement the spectral methods used to solve 
the Navier-Stokes equations. 

Eqs. (II. 1 - 5) are solved by a fractional-step 
spectral method described in detail by Meiron, Orszag 
& Israeli (1981) . The principle changes concern 
the turbulence model, which will now be described. 

The turbulence model used in our calculations is 
that of Jones & Launder (1972) , as modified by Chien 
(1980) . The modeling equations are: 

8U. 8U. 

■'ij = 

11 


D£ _ e 

Dt *^1 ^T k ^9v 


'2 k 


+ 2v- 


-c^y+ 
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'''+ A 
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(II. 7) 
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(II. 8) 
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where the 'eddy viscosity' is given in terms 

of the turbulent kinetic energy k and dissipation 
e by 


V = C 
T V 


-c y 

(1-e ) 


(II. 9) 


here the modeling constants are chosen as 
c^ = 1.35, = c^^d-c^^e 

^21= 1.8, C 22 = 2/9, C 23 = 36, c^ = 0.5 (II. 10) 

c = 0.09, v= 1.57 10 a, = 1, a = 1.3 
V ' k ' e 

2 

where = k /ev is a turbulent Reynolds number 

and the subscript + indicates wall units . The remaining 
modeling constants c^ is chosen a posteriori to 
achieve the best agreement between our calculations 
and available experimental data. 



Results 


In Fig. 1, we plot the skin friction distribution 
as calculated from flow over a flat plate with zero 
gradient. In this figure, \-je. show results obtained 
using our spectral code with both a one equation 
turbulence model and a two-equation model described 
in Sec. II. Results are presented for the two- 
equation model with several choices of the modeling 
constant c^* Chien (1980) suggested the choice 
c^ = 0.0115; we find that c^ = 0.010 gives better 
agreement with experiment. 

The agreement between experiment and numerical 
calculation achieved in Figure 1 is noteable because 
of the apparently significant difference between the 
periodic boundary conditions employed in the computer 
code and the true physical inflow-outflow boundary 
condition. It appears that periodic boundary conditions 
can be used to model the true physics of these boundary 
layer flows with little or no distortion of the results. 
This conclusion is extremely important for time-dependent 
solutions of the Navier-Stokes equations. With inflow- 
outflow conditions, a large computational domain must 
be used to stimulate flow over a long physical body, 
while with periodic boundary conditions, the calculations 
can be done over domain of limited special extent. 


In Figure 2 , we plot the distribution of Reynolds 

2 

stress (normalized by u^ ) as a function of the wall 
variable y^ , again for a flat plate with zero pressure 
gradient. The agreement with the classical data of 
Schubauer is excellent. 

2 

In Figure 3, we plot ev/u^ as a function of 

y_j^, again for a flat plate with zero pressure gradient. 

2 

In Figure 4, we again plot ev/u^ vs y^, at 

21 downstream locations (as determined by the transformation 

X = Ut where U is the free- stream velocity and 

t is the integration time in the quasi-steady state 

achieved by the periodic boundary condition code) . 

These downstream locations have Reynolds numbers 

Rg varying from 2000 to 4400. Observe that the 

results lie on a universal (R. - independent) curve. 

w 

In Figure 5, we plot the kinetic energy distribution 

2 

k (normalized by ) as a function of y_^ for 

flow over a flat plate. The peak of k occurs at 
y_j_ = 25, in good agreement with available measurements. 

We have used the two-equation models together 
with our spectral code to study turbulent flows 
over wavy walls. Here we will report results only 
for flows over sinusoidal wavy walls; 



y = a sin(kx) 


(III.l) 


in the runs reported below, ka = 0.078 so the ratio 
oI amplitude to wave length, a/X, is 0.125. 

In Figure 6,7, we plot-^the variation of 
(downstream velocity in wall units) ''as ^ function 
of y^ in the region of maximum adverse ^nd^^a?vbrable 
pressure gradient, respectively. In Figure 8, we 
plot the distribution of boundary layer thickness 
as a function of position over the wave. In 
agreement with experiment, the bounary layer is 
thickest at the tough of the wave. 

In Figure 9, contours of the turbulent kinetic 
energy distribution k, over the wave surface are 
plotted. In Figure 10, a similar plot of the 
turbulent dissipation , e, is made. From these 
Figures, it is apparent that both turbulent kinetic 
energy and dissipation tend to become smaller in the 
region of favorable pressure gradients. As the 
waves get steeper, this indicates the tendency for 
relaminarization in these regions, so there may be 
a serious problem with two-equation turbulence models 
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in regions of favorable pressure gradient induced be 
steep curvature. 

In Figure 11, we plot contours of the pressure 
distribution over the wavy wall. The maximum value 
of the pressure coefficient Cp and the phase of 
the pressure distribution over the wall 
with available experimental dut-a 


correlations . 



Figure 1 . Variation of coefficient of skin friction for 
a flat plate boundary layer vs R . Here the 
solid curve is obtained from the°tv/o-equation 
model {H.6)“(II,10) with c^ = 0.10. The other 
points are: ^ Model solution with c^ = 0.02; 

V^odel solution with C 3 = 0.0115, as used by 
Chien; □ One equation model; O C 3 = 0.008; 

X Experimental results due to Wexghardt. 
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4 

Figure 4. A plot of ev/u vs at 21 x stations 
with Rg varying from 2000 to 4400. The 
calculations v/ere made for a flat plate 
with zero pressure gradient using the 
two-equation model with c^ = 0.0115. 
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Figure 6. A plot of U_^ vs y_^ for flow over a wavy 

surface. The plot is made at the location 
of maximura favorable pressure gradient with 
ka = 0.0785. 
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CONTOUR FROM 0. TO .29000E-02 

CONTOLfl INTERVAL IS .10000E-D3 LABELS SCALED QT .IOOOOE*06 


Figure 10. A contour plot of the turbulent dissipation for 
flow over a wavy surface with ka = 0.0785. 
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